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Abstract 

The  decomposition  of  unconfined  rigid  polyurethane  foam  has  been  modeled  by  a  kinetic 
bond-breaking  scheme  describing  degradation  of  a  primary  polymer  and  formation  of  a  thermally 
stable  secondary  polymer.  The  bond-breaking  scheme  is  resolved  using  percolation  theory  to 
describe  evolving  polymer  fragments.  The  polymer  fragments  vaporize  according  to  individual 
vapor  pressures.  Kinetic  parameters  for  the  model  were  obtained  from  Thermal  Gravimetric 
Analysis  (TGA).  The  chemical  structure  of  the  foam  was  determined  from  the  preparation  tech¬ 
niques  and  ingredients  used  to  synthesize  the  foam.  Scale-up  effects  were  investigated  by  simu¬ 
lating  the  response  of  an  incident  heat  flux  of  25  W/cm2  on  a  partially  confined  8.8-cm  diameter 
by  15-cm  long  right  circular  cylinder  of  foam  which  contained  an  encapsulated  component.  Pre¬ 
dictions  of  center,  midradial,  and  component  temperatures,  as  well  as  regression  of  the  foam  sur¬ 
face,  were  in  agreement  with  measurements  using  thermocouples  and  X-ray  imaging. 
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EXECUTIVE  SUMMARY 

This  report  describes  the  decomposition  of  unconfined  rigid  polyurethane  foam  using  a 
kinetic  bond-breaking  scheme,  lattice  statistics,  and  vapor  liquid  equilibrium.  Bond-breaking  is 
based  on  degradation  of  a  primary  polymer  and  formation  of  a  thermally  stable  secondary  poly¬ 
mer.  The  bond-breaking  scheme  is  resolved  using  percolation  theory.  The  polymer  fragments 
vaporize  according  to  individual  vapor  pressures.  This  report  describes  the  chemical  structure  of 
the  polyurethane  foam  used  in  this  study,  the  chemical  mechanism  used  to  describe  bond  break¬ 
ing,  the  details  of  the  lattice  statistics  using  percolation  theory,  and  the  vapor-liquid  equilibrium 
model  used  to  determine  the  split  between  liquids  and  gas  formation.  The  report  also  describes 
the  method  used  to  determine  the  kinetic  coefficients,  a  comparison  between  predicted  and  mea¬ 
sured  mass  loss  for  various  TGA  experiments,  an  uncertainty  and  sensitivity  analysis  of  the  PUF 
model,  and  a  comparison  between  the  predicted  shape  of  a  regression  front  and  measured  shape  of 
the  regression  front  using  X-ray  tomography.  The  report  ends  with  a  Summary  and  Conclusions 
section. 

1.  INTRODUCTION 

Rigid  polyurethane  foams  are  used  as  encapsulants  to  isolate  and  support  thermally  sensitive 
components  within  weapon  systems.  When  exposed  to  abnormal  thermal  environments,  such  as 
fire,  various  encapsulated  components  are  designed  to  fail  sequentially.  In  hazards  analysis,  the 
thermal  response  of  the  encapsulated  components  depends  primarily  on  the  behavior  of  the  foam. 
Modeling  foam  decomposition  is  a  difficult  problem  not  only  because  of  the  numerical  challenges 
associated  with  steep  reaction  fronts  but  also  the  difficulty  of  describing  important  chemical  and 
physical  processes,  such  as  non-Newtonian  liquid  flow.  In  the  current  report,  thermal  transport 
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and  decomposition  chemistry  are  discussed  in  detail.  Mass  transport,  species  diffusion,  bubble 
mechanics,  fluid  flow,  and  gravitational  effects  are  beyond  the  scope  of  this  report. 

In  the  past,  simplifying  assumptions  regarding  the  encapsulating  foam  have  been  made  for 
predicting  component  failure.  For  example,  the  decomposition  of  foam  has  been  ignored,  by 
assuming  the  foam  was  not  present  and  adjusting  surface  emissivities  to  match  thermocouple 
data,  or  by  changing  physical  properties  of  the  foam  at  prescribed  temperatures.  Such  approxima¬ 
tions  were  consistent  with  equivalent  approximations  necessary  for  single  processor  calculations. 
With  the  advent  of  massively  parallel  computers,  high-consequence  predictions  of  foam  decom¬ 
position  can  be  made  with  a  more  fundamental  foam  decomposition  model  founded  on  experi¬ 
mental  observations  to  determine  accurate  decomposition  rates,  decomposition  species,  and 
physical  properties  of  the  evolving  solid  residue. 

Computational  models  used  in  hazards  analysis  at  Sandia  National  Laboratories  are 
designed  to  accommodate  mass  loss  associated  with  foam  decomposition.  For  example,  Fig.  1 

shows  a  COYOTE1  finite  element  calculation  of  a  block  of  material  containing  nonreactive  com¬ 
ponents  of  various  shapes  exposed  to  a  constant  energy  flux.  In  this  calculation,  the  encapsulating 
material  was  assumed  to  decompose  by  removing  elements  from  the  computational  domain  based 
on  the  element  exceeding  a  specified  temperature,  without  using  a  realistic  decomposition  mecha¬ 
nism.  Actual  tests  of  polyurethane  foam  exposed  to  abnormal  thermal  environments,  such  as  fire, 
show  the  system  response  to  be  more  complex. 

Several  groups  studying  the  decomposition  of  large  macromolecules  have  employed  sta¬ 
tistical  network  fragmentation  models  to  describe  decomposition  chemistry.  For  example, 

Solomon  and  coworkers  have  implemented  a  computationally  intensive  Monte  Carlo  technique 
to  describe  the  breakup  of  coal.  Grant  et  al.3'5  have  used  pseudo  lattice  structures,  referred  to  as 


10 


Introduction 


Fig.  1  Example  calculation  of  inert  components  encapsulated  in  rigid  polyurethane  foam  ini¬ 

tially  at  100°  C  exposed  to  a  constant  flux  on  the  entire  exposed  surface.  Elements 
were  removed  when  element  temperature  exceeded  150°  C.  Although  foam  regres¬ 
sion  is  shown  as  a  function  of  time,  a  decomposition  model  was  not  used  for  this  cal¬ 
culation.  Figure  used  with  permission  from  Gartling.1 


Bethe  lattices,  to  obtain  closed-form  solutions  of  the  network  statistics  as  derived  by  Fisher  and 

Essam.6  These  closed  form  solutions  parallel  the  determination  of  molecular  weight  distributions 
during  polymer  synthesis  leading  to  the  critical  condition  required  to  form  infinite  polymer  net¬ 
works  referred  to  as  “gels”  by  Flory.7  In  the  present  report,  Flory’s  methods  of  building  polymers 
are  used  to  decompose  polymers  by  assuming  that  closed  rings  or  cycles  cannot  form  and  that 
bond  reactivity  is  independent  of  the  size  of  the  polymer  fragment  containing  the  bond.  Bond 
breaking  is  assumed  to  be  a  random  scission  process  that  describes  the  extent  of  reaction. 
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The  PolyUrethane  Foam  (PUF)  decomposition  model  discussed  in  the  present  report  is 
based  on  three  fundamental  aspects  of  thermal  decomposition:  1)  a  kinetic  bond-breaking  mech¬ 
anism,  2)  lattice  statistics  to  describe  the  evolving  polymer  fragments,  and  3)  vaporization  of  the 

small  polymer  fragments  (oligomers)  with  high  vapor  pressures  as  discussed  by  Fletcher  et  al.4'5 
The  Bethe  lattice  statistical  model  is  referred  to  as  percolation  theory,  since  the  technique  has 
been  used  historically  to  describe  fluid  flow  through  a  network  of  permeable  and  impermeable 
sites.  To  use  percolation  theory  in  a  bond-breaking  mechanism,  the  intact  bonds  correspond  to 
impermeable  sites,  and  the  broken  bonds  correspond  to  permeable  sites.  The  PUF  model  parame¬ 
ters  include  the  Arrhenius  kinetics  -  controlling  bond-breaking,  the  initial  bridge  population,  the 
coordination  number,  and  the  average  site  molecular  weight.  The  initial  bridge  population  is 
assumed  to  be  composed  of  both  strong  bridges  and  weak  bridges  (discussed  further  in  Section  3). 
The  coordination  number  describes  the  connectivity  of  thermally  stable  sites  by  bridges. 

In  the  current  report,  the  term  “polymer”  will  be  used  to  describe  the  infinite  lattice  struc¬ 
ture,  the  term  “polymer  fragments”  refers  to  the  finite  fragments  resulting  from  polymer  decom¬ 
position,  and  the  term  “oligomer”  will  be  used  to  describe  the  finite  fragments  that  evolve  into  the 
gas  phase.  The  term  “bridge”  is  used  to  represent  the  part  of  the  polymer  structure  separating 
thermally  stable  “sites.”  “Bridges”  are  connected  to  “sites”  by  chemical  bonds.  When  a  bond  is 
broken,  a  bridge  is  also  broken.  A  more  detailed  discussion  of  “bridges,”  “sites,”  and  “bonds”  can 
be  found  in  Section  2. 

The  remainder  of  this  report  is  organized  as  follows:  The  next  section  describes  the  chem¬ 
ical  structure  of  the  polyurethane  foam  used  in  this  study.  Section  3  describes  the  chemical  mech¬ 
anism  used  to  describe  bond  breaking.  Section  4  presents  the  details  of  the  lattice  statistics  using 
percolation  theory  as  applied  to  bond  breaking.  Section  5  describes  the  vapor-liquid  equilibrium 
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model  used  to  determine  the  split  between  liquid  formation  and  gas  formation.  Section  6  dis¬ 
cusses  the  method  used  to  determine  the  kinetic  coefficients  for  the  bond-breaking  reactions  as 
well  as  presents  a  comparison  between  predicted  and  measured  mass  loss  for  various  TGA  exper¬ 
iments.  Section  7  presents  an  uncertainty  and  sensitivity  analysis  of  the  PUF  model.  Section  8 
shows  a  comparison  between  the  predicted  shape  of  a  regression  front  and  shape  of  the  regression 
front  using  X-ray  tomography.  Some  general  remarks  close  the  report. 


2.  CHEMICAL  STRUCTURE 

Detailed  information  regarding  the  chemical  structure  of  the  polyurethane  foam  is  required  to  pre¬ 
dict  foam  decomposition  using  lattice  statistics.  The  structural  units  and  resulting  polymeric  net¬ 
work  of  many  synthesized  macromolecules,  such  as  polyurethane  foam,  can  be  inferred  from  the 
starting  materials  and  the  synthesis  of  the  macromolecule.  Confirmation  of  the  structure  is  often 

obtained  using  IR  spectroscopy,  solid-state  NMR,8  and  other  analytical  chemistry  techniques. 

The  most  common  chemical  structural  units  of  the  rigid  polyurethane  foam  and  the  distri¬ 
bution  of  these  structural  units  are  shown  in  Fig.  2  and  were  estimated  using  proprietary  synthesis 
details  and  assuming  equal  reactivity  of  the  hydroxyl  groups.  Various  graphic  representations  of 
specific  chemical  structural  units  are  also  shown  in  Fig.  2.  These  structural  units  have  been 
assembled  to  give  a  representative  picture  of  the  foam.  The  model  foam  shows  the  polyurethane 
structure  as  a  large  matrix  (essentially  infinite)  of  toluene  diisocyanate  groups  connected  by  ali¬ 
phatic  bridges  made  from  trimethylol  propane,  adipic  acid,  diethylene  glycol,  and  small  amounts 
of  phthalic  anhydride. 
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Fig.  2  Most  common  chemical  structural  units  and  hypothetical  chemical  structure  of  rigid 
polyurethane  foam.  The  graphic  symbols  are  composed  of  ingredients  used  to  make 
the  specific  foam. 
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Figure  3  shows  the  three  most  common  structural  units  divided  into  sites  and  bridges.  The 
coordination  number,  a+1,  and  molecular  weights  for  specific  structural  components  are  shown  in 
Fig.  3.  This  information  can  be  used  to  determine  the  average  coordination  number,  as  well  as  the 
average  molecular  weight  of  the  sites  and  bridges  as  discussed  further  in  Section  6.  In  the  present 
report,  the  sites  are  assumed  to  be  primarily  composed  of  trimethylol  propane  units  and  the  dieth¬ 
ylene  glycol  units.  The  bridges  are  composed  of  the  TDI  and  adipic  acid  structures. 

As  the  bonds  break  at  elevated  temperatures,  fragments  of  finite  molecular  weight  are  gen¬ 
erated.  The  fragments,  with  low  molecular  weights  and  corresponding  high  vapor  pressures, 
evaporate  to  form  gaseous  oligomers  as  shown  at  the  bottom  of  Fig.  2.  Higher  molecular  weight 
fragments,  with  low  vapor  pressures,  remain  in  the  condensed  phase. 


3.  KINETIC  MECHANISM 

In  this  study,  the  initial  polymer  structure  begins  to  degrade  between  250  and  350°  C  as  bridges 
between  sites  are  either  broken  or  become  thermally  stable  via  the  evolution  of  light  gases,  similar 

to  the  cellulose  decomposition  model  discussed  by  Grant  et  al.3  and  Shafizadeh.9  At  these  tem¬ 
peratures,  secondary  reactions  between  functional  groups  can  also  form  bridges  between  sites, 
creating  a  secondary  polymer  structure.  The  secondary  polymer  is  assumed  to  decompose  analo¬ 
gously  to  the  primary  polymer,  but  at  temperatures  above  350°  C  and  at  a  different  rate. 

The  PUF  bond  decomposition  model  is  an  extension  of  Grant  et  al.3  bond-breaking 
scheme  and  is  the  first  attempt  to  describe  degradation  of  a  polyurethane  foam  using  percolation 
theory  with  vapor  liquid  equilibrium.  All  of  the  PUF  model  parameters  were  obtained  from  data 
taken  near  ambient  pressures,  and  the  model  is  not  expected  to  perform  well  under  heavy  confine- 
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A)  60%  most  probable  structural  unit  of  a  rigid  polyurethane  foam 


B)  20%  most  probable  structural  unit  of  a  rigid  polyurethane  foam 


C)  10%  most  probable  structural  unit  of  a  rigid  polyurethane  foam 


Fig.  3  Three  most  common  structural  units  of  a  rigid  polyurethane  foam  showing  bound¬ 
aries  used  to  relate  chemical  structure  to  a  Bethe  lattice. 


16 


Kinetic  Mechanism 


ment  where  secondary  reactions  between  the  vapor  and  solid  are  significant.  Currently,  an  exper¬ 
imental  program10  is  being  conducted  to  obtain  additional  information  regarding  the 
decomposition  mechanism  with  attention  focused  on  confinement  issues  and  pressure  effects. 

The  kinetic  scheme,  along  with  a  graphical  description,  is  shown  in  Figs.  4  and  5.  This 
9-step  mechanism  with  11  “species”  describes  four  events:  1)  competition  between  local 
crosslinking  and  side-chain  formation  within  the  primary  polymer,  2)  competition  between 
side-chain  evolution  from  the  primary  polymer  and  the  formation  of  a  secondary  polymer,  3) 
competition  between  local  crosslinking  in  the  secondary  polymer  and  side-chain  formation  within 
the  secondary  polymer,  and  4)  side-chain  evolution  from  the  secondary  polymer. 

The  L,  L,  C-| ,  and  C2  represent  bridges  as  shown  in  Figs.  4  and  5.  A  representative  labile 
bridge  in  the  primary  polymer,  L,  is  the  TDI  structural  component  as  shown  boxed  in  Fig.  4.  The 
L  bridge  is  shown  in  Fig.  4  connecting  two  trimethylol  propane  sites  with  urethane  linkages.  The 
labile  bridge  in  the  secondary  polymer,  L,  is  also  shown  boxed  in  Figs.  4  and  5.  The  L  bridge  is 
composed  of  two  toluene  components  linked  by  a  carbodiimide  group,  ~N=C=N~.  The  L  bridge 
is  connected  to  the  polymer  by  two  urethane  linkages.  Thermally  stable  bridges  that  connect 
sites,  C-|  and  C2,  are  shown  in  Figs.  4  and  5.  The  thermally  stable  bridges  are  actually  bonds  and 
do  not  contain  mass.  Under  certain  conditions,  decomposing  polyurethane  foam  produces  a  char 
or  carbonaceous  residue  that  is  thermally  stable.  In  the  present  report,  the  formation  of  thermally 
stable  bonds  contributes  to  char  formation. 

The  8,  d,  and  gj  represent  side-chains  or  “danglers”  in  the  primary  polymer,  side-chains 
or  “danglers”  in  the  secondary  polymer,  and  various  gas  species,  respectively.  Representative 
“danglers”  are  shown  boxed  in  Figs.  4  and  5.  Reversible  reactions  are  included  to  allow  the  “dan- 
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Fig.  4  Kinetic  mechanism  for  decomposition  of  rigid  polyurethane  foam  with  the  primary 
polymer  decomposition  pathway  highlighted. 
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PUF  bond  breaking  mechanism 
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Fig.  5  Kinetic  mechanism  for  decomposition  of  rigid  polyurethane  foam  with  the  secondary 

polymer  decomposition  pathway  highlighted. 
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glers”  to  reattach  to  the  polymer.  Such  reactions  will  likely  be  more  significant  when  the  foam  is 

degraded  under  confinement.  This  nomenclature  closely  follows  Grant  et  al.35 

The  reaction  sequence  begins  when  a  weak  bond  -  e.g.  the  ether  bond  of  the  carbamate 
group  -  is  broken  to  form  either  a  toluene  diisocyanate  (TDI)  side-chain  (8  in  Fig.  5)  or  evolves  as 
a  gas  (g-|  in  Fig.  4)  with  concurrent  stabilization  of  adjacent  clusters  forming  a  stable  bridge  (C-|  in 
Fig.  4).  The  TDI  side-chains,  referred  to  as  danglers,  may  eventually  evolve  as  light  gas  frag¬ 
ments  (g2  in  Fig.  4)  through  subsequent,  slower  reactions.  A  secondary  polymer  may  evolve  as 
two  TDI  side-chains  react  to  form  a  weak  bridge  (L)  with  the  subsequent  evolution  of  C02 
denoted  as  g3.  The  molecular  weight  of  L  is  twice  the  molecular  weight  of  L  minus  the  molecular 
weight  of  g3.  The  degradation  of  the  secondary  polymer  is  assumed  to  be  similar  to  the  degrada¬ 
tion  of  the  primary  polymer  except  the  labile  bridge  in  the  secondary  polymer  has  a  higher  molec¬ 
ular  weight. 

The  proposed  mechanism  in  Figs.  4  and  5  is  consistent  with  preliminary  experimental 
observation.  For  example,  Fig.  6.A  shows  two  bond  types  that  likely  break  during  thermal 
decomposition:  CD  the  ester  bond  of  the  carbamate  group  and  ®  the  ether  bond  of  the  ester  group. 
Experimental  evidence  of  bond  type  (D  breaking  is  shown  in  Fig.  6.B,  where  an  IR  spectra  of 
decomposition  gases  obtained  using  isothermal  TGA-FTIR  analysis  shows  an  increase  in  the 
~N=C=0  stretch  region  of  the  IR  spectra  (wave  number  of  2250)  from  evolving  gases  sampled  at 
350°  C,  400°  C,  and  450°  C.  Furthermore,  solid-state  NMR  results,  shown  in  Fig.  6.C,  show  the 
depletion  of  the  carbonyls  associated  with  the  adipic  acid  and  isocyanate  structure.  The  NMR 
results  are  consistent  with  the  bonds  labeled  (D  and  ®  breaking.  Although  additional  information 
regarding  decomposition  chemistry  is  evident  in  the  NMR  data,  complete  analysis  was  not  avail- 
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100 


Wave  Number  x  10'3 


Fig.  6  A)  Most  probable  chemical  structure  showing  two  possible  weak  bonds  labeled  CD  and 
©.  B)  FTIR  spectra  showing  evidence  of  urethane  bond  rupture.  C)  Solid-state  NMR 
evidence  showing  change  in  carbonyls  associated  with  adipic  acid  and  isocyanate  struc¬ 
tures.  Other  information  available  in  the  NMR  data  has  not  been  fully  analyzed. 
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able  because  of  insufficient  funding.  The  mechanism  in  Fig.  5  is  also  consistent  with  measured 
decomposition  products,  such  as  TDI  and  trimethoxypropane,  as  determined  from  desorption  tube 

analysis.10 

Table  1  gives  the  mechanism,  rate  equations,  and  initial  conditions  for  the  PUF  model. 
The  general  reaction  rate,  rj,  for  the  PUF  model  bond-breaking  scheme  is  described  by: 

11 

r;  =  ¥r) H”?  >  (!) 

i  =  1 

where  Nt  represents  either  population  parameters  L,  L,  8,  d,  C-|,  C2  or  species  parameters  g-|,  g2, 
g3,  g4,  g5.  The  concentration  matrix  is  represented  by  p;y,  which  is  given  in  the  footnote  of 
Table  1.  The  expressions  for  the  kinetic  coefficients,  kj  (T),  are  given  in  an  Arrhenius  form: 

kj(T)  =  Aj exp(-E j/ RT) ,  (2) 

where  Aj  (1/s),  Ej  (cal/mol  or  J/mol),  and  R  (1.987  cal/mol-K  or  8.314  J/mol-K)  are  the  pre  expo¬ 
nential  factors,  activation  energies,  and  the  universal  gas  constant,  respectively.  The  species  rate 
of  change  is  given  by: 

j 

dN/dt  =  £  V;7r;  ,  i  =  1,  .../,  (3) 

j  =  1 

where  V,y  are  the  stoichiometric  coefficients  of  the  jth  reaction  as  given  in  the  footnote  of  Table  1. 
The  kinetic  mechanism  requires  the  Arrhenius  parameters,  Aj  and  Ej,  to  be  supplied  for  each  reac¬ 
tion.  The  initial  labile  bridge  population  ( L0 )  and  the  initial  stable  bridge  population  (C0)  need  to 
be  specified  as  an  initial  condition.  These  parameters  may  represent  the  extent  of  polymer  curing. 

The  initial  value  problem  described  by  Eqs.  (l)-(3)  is  solved  using  DEBDF11  -  a  variable  order, 
backward  difference,  ordinary  differential  equation  solver  package. 
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Tablel.  Mechanism,  rate  equations,  and  initial  conditions  for  the  PUF  model* 


Rxn 

Mechanism 

Species 

Rate  Equations 

Initial  Conditions 
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L  (0)  =  0 
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*  The  concentration  (jLly)  and  stoichiometric  matrices  (Vy)  with  the  reactions  loaded  into  9  columns  (j  =  9) 
and  the  species  loaded  into  11  rows  (z  =  11)  are  written  as: 
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4.  STATISTICAL  MODEL 

Percolation  theory,  historically  applied  to  describe  fluid  flow  through  a  network  of  permeable 
and  impermeable  sites,  is  used  to  characterize  the  degraded  foam  structure  with  regard  to  the  size 
and  concentration  of  finite  fragments  or  oligomers.  The  population  variables  determined  from  the 
kinetic  mechanism  can  be  used  to  determine  the  fraction  of  sites  that  are  connected  by  bridges. 
For  example,  the  fraction  of  intact  bridges,  p,  can  be  determined  as  follows: 

p  =  JL  +  L  +  Cj  +  C2 .  (4) 

Percolation  theory,  using  Bethe  lattices,  is  computationally  fast,  reproducible,  and  replicates 

results  from  more  versatile,  yet  computationally  demanding,  Monte  Carlo  methods.12  Bethe  lat¬ 
tices  do  not  contain  loop-backs  and  are  similar  to  tree-like  structures.  Monte  Carlo  techniques 
generate  random  numbers  between  0  and  1  for  each  bridge.  The  bridges  are  either  intact  (if  the 
random  number  is  between  0  and  p)  or  broken  (if  the  random  number  is  between  p  and  1).  Monte 
Carlo  techniques  can  be  used  with  realistic  chemical  structures  provided  large  realizations  with 
many  bridges  are  used  to  describe  the  large  macromolecules. 

One  advantage  of  the  Monte  Carlo  technique  over  percolation  theory  is  that  the  polymer 
fragments  can  be  removed  from  the  network  and  quenched  without  affecting  the  decomposition 
kinetics  of  the  parent  foam  structure.  If  the  gas-phase  polymer  fragments,  referred  to  as  oligo¬ 
mers  in  the  present  report,  are  kept  at  the  same  temperature  as  the  parent  foam,  the  oligomers  ther¬ 
mally  crack  to  form  lighter  gases  and  smaller  oligomer  fragments.  The  cracking  phenomena  is 
properly  accounted  for  using  Bethe  lattice  statistics.  In  an  open  system,  oligomers  do  not  decom¬ 
pose  significantly  when  cooled  to  low  temperatures  and  must  be  accounted  for  separately  as  dis¬ 
cussed  by  Kerstein.13  In  the  present  report,  oligomers  are  assumed  to  continue  to  react.  This 
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assumption  does  not  change  the  prediction  of  the  overall  condensed  fraction  which  is  used  as  the 
criteria  for  element  removal  as  discussed  in  Section  7.2. 

A  detailed  formulation  of  percolation  theory  based  on  Bethe  lattices  is  discussed  in  detail 

by  Grant  et  al.3  The  mass  fraction  of  finite  polymer  fragments  produced  from  the  thermally 
degrading  foam  depends  on  the  population  of  intact  bridges,  p,  determined  from  the  kinetic  mech¬ 
anism  and  the  coordination  number,  c  +  1.  The  lattice  statistics  are  independent  of  the  type  of 
bridges  connecting  sites  provided  the  bridges  can  be  distinguished  as  either  broken  or  intact. 
Thus,  percolation  theory  can  be  used  with  any  kinetic  mechanism  provided  that  the  number  of 
broken  or  intact  bridges  is  determined.  A  more  complex  percolation  model  based  on  two  inde¬ 
pendent  networks  has  been  developed  to  allow  reattachment  associated  with  crosslinking.12  The 
“2-c”  model  was  developed  to  allow  the  coordination  number  of  the  pristine  macromolecule  to 

increase  with  extent  of  reaction  due  to  crosslinking.12  In  the  current  work,  the  average  coordina¬ 
tion  number  does  not  change  with  extent  of  reaction,  since  the  sites  are  primarily  composed  of  tri- 
methylol  propane  units.  The  statistics  used  for  the  foam  decomposition  model  in  the  present 
report  are  based  on  a  single  coordination  number. 

Figure  7  shows  Bethe  lattices  with  p  =  1  (all  bridges  intact),  and  p  =  0.67  (67  percent  of 
the  bridges  intact),  for  a  +  1  =  2.2  and  a  +  1  =  4.  The  smaller  coordination  number  (2.2)  is  similar 

to  the  lattice  used  by  Solomon  et  al.14  in  a  more  complex  Monte  Carlo  simulation  used  to  describe 
lattice  statistics.  The  average  coordination  number  of  the  three  most  common  structural  units 
shown  in  Fig.  3  is  2.8. 

As  the  number  of  broken  bridges  increases,  the  fraction  of  finite  polymer  fragments 
increases  relative  to  the  fraction  of  sites  belonging  to  the  infinite  network.  The  infinite  network 

no  longer  exists  below  a  critical  bridge  population  of  1/a  as  shown  by  Fischer  and  Essam.6  The 
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(a  +  l)ave  =  2.2 
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Fig.  7  Bethe  lattices  with  coordination  numbers  2.2  and  4.0.  Each  coordination  number  ex¬ 
ample  is  shown  with  critical  bridge  populations  of  1.0  and  0.67. 


critical  bridge  populations  for  the  two  Bethe  lattices  shown  in  Fig.  7  are  0.83  and  0.33  for  coordi¬ 
nation  numbers  of  2.2  and  4,  respectively.  Thus,  the  a  +  1  =  2.2  lattice,  with  p  =  0.67,  is  com¬ 
posed  of  finite  polymer  fragments  exclusively;  the  higher  coordination  lattice  with  p  =  0.67  still 
maintains  the  infinite  polymer  network. 

For  Bethe  lattices,  the  number  of  finite  polymer  fragments  can  be  determined  from  the 
coordination  number  and  the  bridge  population,  p.  The  probability,  Fn,  that  any  given  site  is  a 
given  member  of  finite  polymer  fragment  of  n  sites  with  s  bridges  is 


Fn  =  anp\\-p)\ 

(5) 

where 

s  =  n  -  1; 

(6) 

and 

x  =  n(a  -  1)  +  2. 

(7) 
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x  is  the  number  of  broken  bridges  on  the  perimeter  of  the  polymer  fragment  with  ^-bridges  con¬ 
necting  n-sites  as  shown  for  the  trimer  in  Fig.  8  The  broken  bridges  isolate  the  finite  polymer 
fragments  from  the  infinite  lattice.  Figure  9  shows  the  weight  fraction  of  various  finite  polymer 
fragments  as  a  function  of  the  bridge  population  for  a  =  2.  The  number  of  different  ways  to  form 

such  fragments  is  represented  by  an :3 


f  <3  +  1 


r(nct  +  2) 


Uo+ur(n)r[n(o-i)  +  3] 


(8) 


where  T  represents  the  standard  gamma  function. 

Equation  (5)  describes  the  probability  that  any  given  bridge  belongs  to  an  n-mer.  The 
probability  of  belonging  to  an  n-mer  is  the  probability  that  the  given  bridge  is  intact  or  occupied 

( ps ),  multiplied  by  the  probability  that  the  nearest  neighbor  bridges  are  broken  or  unoccupied 
( 1  -p)x.  an  accounts  for  the  distinct  number  of  configurations  possible  for  the  n-mer.  Equations 
(5)-(8)  are  discussed  in  more  detail  in  References  3,  6,  and  7. 

The  total  fraction  contained  in  all  of  the  finite  clusters,  F,  is3 


F=  £  F„  = 

n-  1 


where  p*  is  the  root  of  the  following  equation: 


Vl-pl 

a  +  1 

~p*~ 

l  _  p*_ 

-  p  - 

(o+l)/(o-l) 


z  =  1  =  p(l-p)°  !,or 


p*(l-p*)G  1  -  p(l  -  p)°  1  =  0; 


(9) 


(10) 

(11) 


and  p  is  determined  from  Eq.  (4).  The  fraction  of  sites  contained  in  the  infinite  cluster,  /,  is 


/  =  1  -F. 

labeled  in  Fig.  9  as  the  infinite  cluster,  °o-mer. 


(12) 
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Trimer  (oligomer  with  3  sites) 
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Fig.  8  Trimer  showing  a,  n,  s,  and  x. 
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Fig.  9  Weight  fraction  of  finite  polymer  fragments  and  the  infinite  cluster  vs.  bond  population 
for  cr  =  2.  The  curves  have  been  calculated  from  Eqs.  (5)  and  (12). 
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For  every  permissible  value  of  z  in  Eq.  (10),  two  roots  exist  in  the  range  0<p<  1 ,  except  at 
the  maximum  where  the  roots  are  identical  as  can  be  seen  in  Fig.  10.  When  G  =  2,  z  is  a  simple 
quadratic  function  with  a  maximum  at  0.5.  A  critical  bridge  population  defining  the  existence  of 
an  infinite  lattice  can  be  determined  by  setting  the  first  derivative  of  z  to  0  and  solving  for  p. 

^  =  0  =  (1  -ap)(l  -p)°-2  =>  p  =  I  (13) 

dp  a 

For  a  cluster  to  propagate  indefinitely  from  a  given  intact  bridge,  the  adjoining  site  must  contain 
at  least  one  other  intact  bridge.  Since  p  is  the  probability  that  a  bridge  is  intact  and  G  bridges  radi¬ 
ate  from  the  adjoining  site,  an  infinite  cluster  exists  only  when  Gpc  >  1  or  pc>  1/g.  The  trivial 
root  of  Eq.  (1 1),  p*  =  p,  should  be  used  when  p  is  less  than  pc,  making  F  =  1.  The  other  root  that 
is  bounded  by  0 <p*<pc  should  be  used  when  p  is  greater  than  pc. 


G  =  2  G  =  3 


Fig.  10  Plot  of  Eq.  (11)  for  various  G,  p*,  and  p. 
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The  mass  fraction  and  molecular  weight  of  each  polymer  fragment  bin  can  be  determined 
by  relating  the  total  mass  and  mass  associated  with  finite  polymer  fragments  on  a  site  basis.  A 
site  is  defined  as  the  portion  of  the  polymer  that  can  be  isolated  by  bridges  as  shown  enclosed  in 
dashed  lines  in  Fig.  1 1 .  The  total  mass  per  site  is: 

=  l-c„)(o+l),  (14) 

where  Ma  is  the  molecular  weight  of  the  site  as  represented  by  a  sphere  (or  mer)  as  shown  in 
Fig.  11,  Mb/2  is  half  the  molecular  weight  of  the  bridges  connecting  the  spheres,  and  ( 1  -ca)  repre¬ 
sents  the  fraction  of  bridges  which  are  labile,  and  G+l  is  the  coordination  number. 

The  mass  of  gas  released  can  also  be  expressed  on  a  site  basis  as 


mg  = 


vfM 
2 


b\  (  b  —  ^  g  3 

-T  <9, +  92)  +  '— - s- 


2  ,(g4  +  g5)  +  (~2  )(93> 


(G+l), 


(15) 


where  Mb  is  the  molecular  weight  of  gj  and  g2,  2  Mb  -  Mg3  is  the  molecular  weight  of  g4  and  g5, 
and  Mg3  is  the  molecular  weight  of  g2.  These  molecular  weights  are  divided  by  two  to  normalize 
the  gas  populations  to  the  total  mass  per  site  given  in  Eq.  (14)  by  noting 


“charred” 
bridge 


labile 

bridge 


1 

side-chain 
dangler” 


Fig.  1 1  5-mer  showing  sites,  bonds,  and  danglers. 
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9i  +  92  +  03 +  204  +  2g5  ->(l-c0)  as  time— (16) 
The  mass  of  a  finite  polymer  fragment  expressed  on  a  site  basis,  mn,  can  be  determined 

from  the  mass  of  the  finite  fragment,  Mn,  multiplied  by  the  n-site  polymer  fragment  population  on 

a  site  basis,  Qn : 


(17) 

(18) 

(19) 

Mn  is  the  molecular  weight  of  the  «,/!-polymer  fragment  bin.  Equation  (18)  and  other  subsequent 
boxed  equations  are  primary  equations  used  in  the  PUF  model.  The  first  term  in  Eq.  (18)  repre¬ 
sents  the  number  of  mers  in  the  n-mer  multiplied  by  the  mer  molecular  weight.  The  second  term 
in  Eq.  (18)  represents  the  number  of  bridges,  n  -  1,  in  the  polymer  fragment  multiplied  by  the 
mass  of  the  bridges  of  type  L.  The  third  term  in  Eq.  (18)  represents  the  number  of  bridges  in  the 
polymer  fragment  multiplied  by  the  mass  of  the  bridges  of  type  L.  The  last  two  terms  in  Eq.  (18) 
represent  the  weight  of  the  side-chains,  or  “danglers,”  which  can  evolve  over  time.  The  factors 
8/(1  -p)  and  d/(l -p)  represent  the  fraction  of  side-chains  of  type  8  and  d,  respectively.  The  frac¬ 
tion  of  broken  bridges  with  one  side-chain  being  formed  from  each  broken  bridge  is  represented 
by  (1  -p),  and  x  is  the  number  of  bridges  that  isolate  a  polymer  fragment  as  given  in  Eq.  (6).  The 
factor  “2”  in  the  denominator  of  the  last  two  terms  in  Eq  (18)  is  needed  to  be  consistent  with  Eqs. 
(14)  and  (15). 


mn  =  MnQn 


where, 


M„  =  nMa  +  <«  -  1  )M„ |  £)  +  (n  -  1  )(2 M„  -  Mg3)^  + 


Mi 


t8 


.(1  -p). 


+ 


(2  Mb-M„,) 


xd 


(1  -p). 


a,  =  Fjn  =  by-'a-pf*-"*2. 
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The  mass  of  all  finite  polymer  fragments,  my,  can  be  determined  by  summing  all  of  the 
polymer  fragments  as  follows: 


mf  =  £  mn  =  £  MnQn  •  (20) 

n  -  1  n  -  1 

The  infinite  summation  in  Eq.  (19)  can  be  evaluated  by  using  the  following  functions  derived  by 
Grant  et  al.:3 


K  =  I  Q,  = 

n-  1 


l  - 


a  +  1 


P__ 
.  P 


*■ 


(o+l)/(o-l) 


(21) 


F  =  £  nQn  = 

n  =  1 


-P  - 


(22) 


Note  that  F  was  defined  previously  in  Eq.  (9).  Combining  Eqs.  (6),  (18),  and  (20)-(22)  gives 


my  =  MaF<P  +  MbK£l , 


where 


.  ,  L\  r8(a-l)  RL  Rd(a-1) 

<P  =  1  +  r\  -  +  — h. - {  +  —  +  — ) - ^ , 

PJ  2(1  -p)  p  2(1  -p) 


„  8  (L\  R'd  R'L 

Q.  =  — - -  -  -  Gl¬ 


and 


(1  -p)  \p)  (1  -p)  P 

r  =  Mb/Ma, 

R  =  (2  MB-Mg3)/MA, 

R'  =  (2  MB-Mg3)/MB. 


(23) 

(24) 

(25) 

(26) 

(27) 

(28) 


The  mass  fraction  of  gas,  finite  polymer  fragments,  and  infinite  lattice  can  be  calculated 
from  Eqs.  (14),  (15),  and  (23)  as  follows: 


fe  =  ™e/mt  = 


[r(#i  +  8 2 )  +  R  83  +  F(gA  +  g5)](cr  +1) 
2  +  r(l-c0)(o+l) 


(29) 
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ft  =  mf/mt  =  .  tz  2  ..  —[<PF  +  £lrK]  (30) 

foo  =  1-fg-ff,  (31) 

where  the  infinite  lattice  fraction  is  referred  to  as  the  “gel  fraction”  by  Flory.7  In  Eq.  (29),  r  and  R 
were  defined  previously  in  Eqs.  (26)  and  (27),  respectively.  The  variable,  R",  represents  the  ratio 
of  the  molecular  weight  of  g3  to  the  molecular  weight  of  a  site,  R"  =  Mg3/MA.  The  mass  fraction 

of  the  finite  oligomers,  fp  is  the  same  as  given  by  Grant  et  al.3  except  for  the  definitions  of 
<E>  and  Q.. 

The  primary  variables  of  interest  for  the  network  statistical  model  are  molecular  weight  of 
the  various  gases  (Mb  is  the  molecular  weight  of  g  |  and  g2;  2  Mb  -  Mg3  is  the  molecular  weight  of 
g4  and  g5;  and  Mg3  is  the  molecular  weight  of  g3),  the  mass  fraction  of  the  gas  (Eq.  29),  the 

molecular  weight  of  the  n,/!-polymer  fragment  (Eq.  18),  and  the  mass  fraction  of  the  n^-polymer 
fragment.  Equation  (30)  gives  the  mass  fraction  of  all  finite  fragments.  The  mass  fraction  for  the 

individual  nth- fragments  can  be  determined  from  Eq.  (30)  by  replacing  F  and  K  by  nQn  and  Qn, 
respectively: 


fn  = 


2  +  r(l  - c  )(o  +  1) 


[<D  nQn  +  SlrQn]. 


(32) 


A  separate  model  must  be  used  to  determine  the  amount  of  finite  polymer  fragments  trans¬ 
ported  to  the  gas-phase.  The  oligomers  that  are  transported  into  the  gas-phase  are  sometimes 

referred  to  as  tar.3,15  Evolution  of  tar  involves  condensed-phase  diffusion  to  the  surface, 
vapor-liquid  equilibrium  at  the  surface,  and  diffusion  through  a  gaseous  boundary  layer.  The 
finite  polymer  fragments  that  remain  in  the  condensed  phase  are  sometimes  referred  to  as  meta- 

plast.15  The  metaplast  (finite  polymer  fragments)  and  gel  (infinite  polymer)  are  referred  to  as 
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char.15  Ideal  vapor- liquid  equilibrium  is  used  as  a  first  order  approximation  to  the  complex  tar 
diffusion  process. 


5.  VAPOR-LIQUID  EQUILIBRIUM  MODEL 

The  fraction  of  the  gaseous  oligomers  can  be  determined  by  using  a  simple  vapor-liquid  equilib¬ 
rium  relationship  using  a  combination  of  Dalton’s  law  and  Raoult’s  law.  In  the  current  report,  a 
standard  multicomponent  isothermal  flash  calculation  was  used  to  determine  the  split  between 

vapor  and  condensed  phases  following  the  procedure  used  by  Fletcher  et  al.5 

In  the  nineteenth  century,  Dalton  observed  that  the  pressure  exerted  by  a  mixture  of  ideal 
gases  is  the  same  as  the  sum  of  pressure  exerted  by  the  individual  gases  occupying  the  same  vol¬ 
ume.16  Dalton’s  law  states  that  the  partial  pressure  is  the  gas-phase  mole  fraction,  yt,  multiplied 
by  the  total  pressure,  P: 

Pi  =  JiP-  (33) 

Raoult  observed  that  the  partial  pressure  in  the  vapor,  Pt,  of  each  component  in  an  ideal 

liquid  solution  was  proportional  to  the  components  mole  fraction  in  the  liquid,  x(,  with  the  propor¬ 
tionality  constant  being  the  vapor  pressure  of  the  pure  component  at  the  system  pressure,  P*:17 

Pi  =  Pi**i-  (34) 

Ideal  solutions  with  chemically  similar  components  obey  Raoult’s  law,  and  dissimilar  species 
deviate  strongly  from  Raoult’s  law.  Eqs.  (33)  and  (34)  can  be  combined  to  define  a  “K- value”  or 

vapor-liquid  equilibrium  ratio:17 
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Ki  =  y->  =  p-. 

1  X ■  P 


(35) 


1 1 

A  standard  multicomponent  isothermal  flash  calculation  can  be  used  with  Eq.  (35)  to 
determine  the  split  between  vapor  and  liquid.  A  mole  balance  between  the  initial  condensed 
phase,  F,  and  the  resulting  vapor  phase,  V,  and  liquid  phase,  L,  is 


F  =  V  +  L,  (36) 

where  F,  V,  and  L  refer  to  total  moles  of  feed,  vapor,  and  liquid,  respectively.  The  feed  in  the 

PUF  decomposition  model  includes  the  sum  of  the  finite  polymer  fractions  determined  from  per¬ 
colation  theory  on  a  mole  basis  determined  from  Eqs.  (18)  and  (32).  The  feed  also  includes  the 
gases  gj  through  g5  that  are  calculated  with  the  kinetic  mechanism.  The  moles  of  the  individual 
components  in  the  feed,/),  vapor,  vy,  and  liquid  /,-,  sum  to  F,  V,  and  L,  respectively: 


n  n  n 

F  =  E/f- v  =  £v,.,  and  L  =  £/,  (37) 

i  =  1  i  =  1  i  =  1 

fi  =  Z1F ’  vi  =  yiv ’  and  h  =  xiL-  (38) 

Z(,  yi,  and  xt  represent  the  mole  fraction  of  the  ith  component  in  the  feed,  vapor,  and  liquid,  respec¬ 
tively.  An  expression  for  .iy  can  be  determined  from  a  single  component  material  balance: 


fi  =  vi  +  h 

ztF  =  KjXjV  +  xtL  =  KiXiV  +  x.[F-V]  =  xfF  -V  +  KiV), 


or  = 


1  (F-V  +  K-V) 


i-V  +  K‘v 


F  F 


1)^+1 


An  expression  of  >y-  can  be  determined  from  Eqs.  (35)  and  (41)  as 


(39) 

(40) 

(41) 
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»■  =  Kixi  = 


K*i 


(42) 


A  stable  implicit  equation  to  solve  for  V7F  can  be  determined  using  the  identity  L  Xj  =  L  yj  =  1,  or 


n  n 


E  »  -  E  *,•  =  °  =  E 


;  =  l  /  =  l 


y  ^  y 

/=  l(ATf-  1)^;+  1  1=1  (Kj  —  1 )  ~  +  1 


(43) 


or 


0=  E 


y 


(44) 


/=i(^.-l)^  +  l 


V/F  is  determined  iteratively  using  the  zeroin  solver.18  Eq.  (44)  is  the  Rachford-Rice  equation.17 

The  gases  gj  through  g5  must  be  combined  with  the  finite  fragments  to  calculate  the  feed 
mole  fractions,  z(.  The  “K- values”  defined  in  Eq.  (35)  can  be  determined  from  the  vapor  pressure 


of  the  pure  components,  P*,  divided  by  the  system  pressure,  P.  Fletcher  et  al.5  have  proposed  the 
following  vapor  pressure  correlation  for  high  molecular  weight  organic  molecules: 


Pn*  =  87100exp 


r  -299M 

T 


0.590\ 

n 


(45) 


V  J 

where  Pn*  is  the  vapor  pressure  of  the  pure  n^-polymer  fragment  in  atmospheres,  Mn  is  the 

molecular  weight  of  the  n^-polymer  fragment  in  g/mol  as  determined  from  Eq.  (18),  and  T  is  the 
temperature  in  K.  The  coefficients  in  Eq.  (45)  were  obtained  from  coal  tars  with  molecular 
weights  ranging  from  110  to  315  amu.  The  vapor  pressure  predicted  with  Eq.  (45)  increases  with 
lower  molecular  weight  species.  The  functional  form  of  Eq.  (45)  is  similar  to  the  Clausius-Clap- 
eyron  equation  which  is  usually  used  to  determine  the  vapor  pressure  at  any  temperature  when  the 
vapor  pressure  is  known  at  another  temperature: 
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P*  =  P*knownQX  P 


-A 


R  \T  T 


known 


(46) 


where  A Hv  is  the  molar  heat  of  vaporization  and  R  is  the  gas  constant. 

Fletcher  et  al.5  have  shown  acceptable  agreement  with  boiling  point  data  for  111  organic 
compounds  at  pressures  of  0.007,  0.08,  1,  and  10  atm  using  the  correlation  given  in  Eq.  (45).  The 
111  organic  compounds,  with  molecular  weights  as  high  as  244,  did  not  contain  long  chain 
alkanes  or  compounds  with  more  than  two  oxygen  atoms.  Long  chain  hydrocarbons  are  not 
expected  to  occur  in  the  polyurethane  tar  in  significant  quantities.  However,  some  of  the  expected 

products  measured  by  Chu  et  al.,10  such  as  trimethoxypropane,  do  contain  more  than  two  oxygen 
atoms. 


The  vapor  pressures  of  expected  polyurethane  decomposition  products,  such  as  alcohols, 
acids,  esters,  and  isocyanates,  are  likely  to  behave  according  to  the  correlation  in  Eq.  (45)  within 
each  family  of  products.  The  most  prevalent  decomposition  products  are  TDI  (see  g-|  and  Q2  in 
Fig.  4)  and  carbon  dioxide  (C02).  In  the  current  work,  specific  vapor  pressure  correlations  for 
TDI  and  C02  were  used  rather  than  the  correlation  in  Eq.  (46).  Specifically,  the  TDI  vapor  pres¬ 
sure  was  used  for  g-|  and  g2;  the  C02  vapor  pressure  was  used  for  g3;  and  the  correlation  given  in 
Eq.  (45)  was  used  for  g4,  g5,  and  the  polymer  fragments  determined  from  percolation  theory. 
More  experiments  are  needed  to  improve  the  vapor  pressure  correlation  for  polyurethane  foam 
fragments. 


6.  ESTIMATING  THE  PUF  DECOMPOSITION  MODEL  PARAMETERS 


Nominal  values  of  the  parameters  used  in  the  PUF  model  are  shown  in  Table  2  with  an  estimate  of 
the  high  and  low  values  used  in  the  sensitivity  analysis  discussed  in  Section  7.1.  The  initial  bridge 
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population  parameters,  pa  and  ca,  were  estimated  based  on  the  degree  of  polymer  curing.  The 
coordination  number  and  molecular  weights  were  calculated  from  the  three  most  common  struc¬ 
tural  units  of  the  foam  as  shown  in  Figs.  2  and  3.  Figure  3  shows  sites  made  from  the  trimethylol 

propane  units  and  the  diethylene  glycol  units.  The  bridges  are  composed  of  TDI  and  adipic  acid 

structures.  Table  3  shows  how  the  average  molecular  weights  and  coordination  numbers  for  the 
60%,  20%,  and  10%  most  probable  structural  units  were  calculated.  The  molecular  weights  and 
coordination  number  in  Table  2  were  determined  by  taking  the  probability  normalized  molecular 
weights  and  coordination  numbers  from  Table  3  as  follows: 

Ma  =  s|(131>  +  |i(118>  +  ^(104'  =  125g/mol  (47) 

Mb  =  jj|(163)  +  2|(160)  +  2i(155)  =  161  g/mol  (48) 

0+1  =  sl(3)+li(2'5)+li<2)  =  28 '  <49) 

The  Arrhenius  parameters  were  obtained  by  minimizing  the  root  mean  squared  error 
between  the  calculated  and  measured  mass  loss  for  six,  unconfined,  isothermal  TGA  experiments 
using  DAKOTA  (Design  Analysis  Kit  for  OpTimizAtion,  Ref.  20).  The  unconfined  TGA  experi¬ 
ments  considered  nominally  5-mg  samples  in  open  platinum  pans  and  were  designed  to  examine 
decomposition  mechanisms  under  conditions  that  minimize  mass  transfer  and  reversible  and  sec¬ 
ondary  reactions.10,21  Figures  12.A  and  12.B  show  a  comparison  between  predicted  and  mea¬ 
sured  condensed  mass  fractions  (foam  mass  measured  by  the  TGA  divided  by  the  initial  foam 
mass)  for  the  isothermal  TGA  experiments  used  to  obtain  the  kinetic  parameters.  The  single  iso¬ 
thermal  TGA  samples  were  quickly  ramped  (-10-15  min)  to  the  temperatures  indicated  in  Fig. 
12. A.  The  dual  isothermal  TGA  samples  were  ramped  to  300°  C  and  held  for  about  two  hours, 
then  ramped  to  the  temperatures  indicated  in  Fig.  12.B.  The  three  single  isothermal  experiments 
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Table  2.  PUF  model  parameters  with  estimated  values  for  sensitivity  analysis* 


Variable 

Symbol 

Description 

Nominal 

Value 

Low  Value  (+) 

High  Value  (-) 

Po 

Initial  bond  population^ 

0.85 

0.8 

0.9 

co 

Initial  strong  bond  population^ 

0.10 

0.05 

0.15 

G+l 

Coordination  number^ 

2.8 

2.1 

3.0 

Mc 

Site  molecular  weight^ 

125  g/mol 

104  g/mol 

131  g/mol 

Mb 

Bridge  molecular  weight^ 

161  g/mol 

155  g/mol 

163  g/mol 

Mgs 

Molecular  weight  of  g3^ 

44  g/mol 

44  g/mol 

44  g/mol 

p 

Pressure^ 

83  KPa 

82  KPa 

84  KPa 

El 

Activation  Energy  for  reaction  1 : 

53400.  cal/mol 

0.9  xE2 

l.lxEj 

e2 

Activation  Energy  for  reaction  2^ 

44580.  cal/mol 

0.9  xE2 

1.1  xE2 

E3 

Activation  Energy  for  reaction  3*' 

42520.  cal/mol 

0.9  xEs 

l.lxEj 

e4 

Activation  Energy  for  reaction  A 

45900.  cal/mol 

0.9  xE4 

\AxE4 

e5 

Activation  Energy  for  reaction  5f 

44600.  cal/mol 

0.9  xE5 

1.1  xE5 

e6 

Activation  Energy  for  reaction  6  ' 

58060.  cal/mol 

0.9  xE6 

1.\xE6 

E? 

Activation  Energy  for  reaction  7  ' 

51850.  cal/mol 

0.9  XEy 

1.1  XEy 

e8 

Activation  Energy  for  reaction  8  ' 

52830.  cal/mol 

0.9  xE8 

1AxE8 

Eg 

Activation  Energy  for  reaction  9  ' 

56660.  cal/mol 

0.9xEg 

lAxEg 

*  All  frequency  factors,  A,  were  assumed  to  be  equal  to  3.0xl015  s"1. 
f  High  and  low  values  guessed. 

$  High  and  low  values  from  the  10%  and  60%  most  probable  chemical  structure. 


Table  3.  Molecular  weight  and  coordination  number  of  common  foam  structural  units 


Structure 

Ma,  g/moll 

Mb,  g/mol 

a  + 1 

60%  most  probable 

131 

[(4x176) +  112]/5  =  163 

3 

20%  most  probable 

(104+131)/2  =118 

[(3x176)  +  112]/4=  160 

(2  +  3)/2  =  2.5 

10%  most  probable 

104 

[(2x176)  +  112]/3  =  155 

2 

*  See  Fig.  3  for  the  boundaries  used  to  decompose  the  most  probable  structural  units  into  bridges  used 
to  determine  MB  and  sites  used  to  determine  MA  as  well  as  a  +  1 . 
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0  50  100  150  200  250  300  350 


Time,  min 

Fig.  12  Comparison  of  PUF  predicted  (solid  line)  and  measured  (dashed  line)  solid  mass  fraction 
for  A)  single  isothermal  TGA,  B)  dual  isothermal  TGA,  and  C)  three  nonisothermal  TGA 
experiments. 


40 


Estimating  the  PUF  Decomposition  Model  Parameters 


shown  in  Fig.  12.  A  were  intended  to  isolate  the  kinetics  for  the  decomposition  of  the  initial  poly- 

mer  and  formation  of  the  secondary  polymer  as  discussed  in  more  detail  by  Erickson  et  al.  The 
final  three  isothermal  experiments  shown  in  Fig.  12.B  were  ramped  to  300°  C  and  held  until  mass 
loss  was  minimal  and  then  ramped  to  the  final  isothermal  temperature  of  380°  C,  400°  C,  and  420° 
C,  respectively.  The  intent  of  these  final  three  experiments  was  to  isolate  the  kinetics  of  the 
decomposition  of  the  secondary  polymer. 

The  predicted  solid  fraction  is  greater  than  the  measured  solid  fraction  for  the  270°  C  iso¬ 
thermal  TGA  results  shown  in  Fig.  12. A.  The  difference  between  the  predicted  and  measured 
solid  fraction  is  not  as  pronounced  at  300°  C  and  330°  C.  At  the  lower  temperature,  alternate  reac¬ 
tion  pathways  involving  evolution  of  a  site,  such  as  a  trimethylol  propane  with  subsequent 
crosslinking  of  the  resulting  urethane  danglers,  may  explain  the  differences.  Future  work  will 
address  this  issue. 

The  kinetic  parameters,  determined  from  the  six  isothermal  TGA  experiments,  are  given 
in  Table  2.  Three  nonisothermal  experiments,  not  used  to  estimate  kinetic  parameters,  were  simu¬ 
lated  to  validate  the  selection  of  the  kinetic  parameters  for  the  PUF  model.  Fig.  12.C  shows  good 
agreement  between  predicted  and  measured  condensed  mass  fraction  for  the  three  nonisothermal 
TGA  experiments  with  heating  rates  of  5,  20,  and  50  °C/min. 

Figure  13  shows  the  predicted  solid  fraction  and  measured  solid  fraction  for  the  20  °C/min 
TGA  case  for  a  4.687-mg  sample  and  a  13.765-mg  sample.  Differences  in  the  experimental 
results  are  probably  related  to  mass  transport  effects.  However,  both  sets  of  data  show  an  initial 
rapid  decomposition  followed  by  a  slower  rate  of  solid  mass  loss. 

Figure  14  shows  solid  mass  fraction,  population  parameters,  the  gas  and  oligomer  mass 
fractions,  and  gas  molecular  weight  for  the  5  °C/min  ramped  TGA  experiment.  Between  0  and  40 
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Fig.  13  Comparison  of  PUF  predicted  (solid  line)  and  measured  (dashed  and  dotted  lines) 
solid  mass  fraction  for  the  20  °C/min  TGA  experiment.  The  dashed  lines  are  for  a 
4.687-mg  sample,  and  the  dotted  line  is  for  a  13.765-mg  sample. 


minutes  in  Fig.  14.B,  P  is  shown  to  be  the  sum  of  L  and  c-|.  At  40  minutes,  8  starts  to  decrease 
with  a  corresponding  increase  in  L,  a  result  of  side-chain  cross-linking  of  8  to  form  the  bridge  L. 
At  50  minutes,  the  secondary  bridge,  L,  begins  to  form  as  a  result  of  reaction  between  two  8 
side-chains.  The  increase  in  both  L  and  L  causes  the  bridge  population,  p,  to  increase  between  40 
and  50  minutes.  Decay  of  the  weak  bridge,  L,  is  shown  as  the  side-chain  population  8  increases. 
The  strong  bridge  formation  in  the  primary  and  secondary  polymer  are  minimal.  However,  the 
formation  of  strong  bridges  may  become  significant  at  elevated  pressures. 
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Fig.  14  A)  Comparison  of  PUF  predicted  (solid  line)  and  measured  (dashed  line)  solid  mass 
fraction  for  a  5  °C/min  TGA  experiment.  B)  Predicted  population  parameters  for  the 
5  °C/min  TGA  experiment  given  in  A.  C)  Predicted  average  volatile  molecular  weight 
and  various  oligomer  mass  percents  for  the  5  °C/min  TGA  experiment  given  in  A. 
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7.  UNCERTAINTY  AND  SENSITIVITY  ANALYSIS  OF  THE  PUF  MODEL 

Two  techniques  were  used  to  analyze  parameters  in  the  decomposition  model:  1)  a  variance  anal¬ 
ysis  discussed  in  Section  7.1  and  2)  a  response  derivative  analysis  discussed  in  Section7.3.  A  grid 
refinement  analysis  is  also  reported  in  Section  7.2  to  establish  that  solutions  are  not  dependent  on 
the  size  of  elements  used  in  the  simulations. 

The  variance  analysis  is  used  to  examine  primary  PUF  model  parameters  that  effect  the 
prediction  of  condensed  mass  fraction  for  the  20  °C/min  ramped  TGA  experiment.  Parameters 
with  primary  effects  were  the  initial  bridge  population,  P0;  the  lattice  coordination  number,  G+l; 
molecular  weight  of  the  bridging  structures,  Mwb\  and  several  activation  energies  associated  with 
the  formation  of  a  secondary  polymer,  £). 

Two  response  derivative  analysis  are  performed  for  fast  (0.97±0.14  cm/min)  and  slow 
(0.26±0.05  cm/min)  recession  rate  conditions.  For  the  fast  recession  rate  analysis,  the  most  sensi¬ 
tive  parameters  are  the  activation  energies  associated  with  the  formation  of  the  secondary  poly¬ 
mer.  For  the  slow  recession  rate  analysis,  the  most  sensitive  parameters  are  associated  with  the 
thermophysical  properties  of  the  foam  and  the  destruction  of  the  secondary  polymer. 


7.1  VARIANCE  ANALYSIS 

To  determine  the  primary  effects  of  various  PUF  model  parameters,  a  simple  sensitivity  analysis 

is  performed  on  15  parameters  using  a  two-level  20-run  Plackett  Burman22  analysis.  Typically, 
Plackett  Burman  analyses  are  used  in  experimental  design  techniques.  These  techniques  are  also 
useful  in  determining  parameter  sensitivity.  Primary  effects  are  determined  by  comparing  calcu- 
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lated  factor  effects  with  probability  points  of  a  2-sided  t-distribution  with  4  degrees  of  freedom. 
One  degree  of  freedom  is  used  to  estimate  the  mean.  The  primary  response  variable  for  this  study 
is  chosen  to  be  the  condensed  fraction  remaining  during  decomposition  of  polyurethane  foam 
ramped  at  20  °C/min  as  in  the  thermogravimetric  analysis  (TGA).  The  15  parameters  used  for 
this  study  are  given  in  Table  2.  The  sensitivity  of  Mg3  was  not  determined,  since  g3  is  assumed  to 
be  C02.  The  20-run  Plackett  Burman  design  used  in  this  study  is  given  in  Ref.  22. 

To  determine  significant  factor  effects  using  Plackett  Burman  analyses,  a  low  value  and 
high  value  for  each  of  the  PUF  variables  are  selected  as  shown  in  Table  2.  Some  of  the  values  are 
arbitrarily  chosen.  For  example,  the  high  and  low  values  for  the  activation  energies  are  taken  as 
the  nominal  value  ±10%.  The  molecular  weight  values  are  taken  from  the  10%  and  60%  most 
probable  structures.  Twenty  simulations  of  the  TGA  experiment  with  a  20  °C/min.  ramp  are  run 
with  the  parameters  set  to  either  the  low  value  or  the  high  value  as  described  in  Ref.  22.  The  con¬ 
densed  mass  fraction  at  various  times  is  used  as  the  response  variable.  From  the  response  at  these 
times,  a  factor  effect  is  calculated  using  methods  discussed  in  Ref.  22.  The  significant  factors 
based  on  a  90%  confidence  interval  are  shaded  and  boxed  in  Table  4.  The  factor  effects  before 
480  seconds  are  zero,  since  the  temperature  is  below  significant  reaction  thresholds.  Additional 
information  regarding  statistical  inference  using  Plackett  Burman  analysis  can  be  found  in  Refer¬ 
ence  22. 

The  Plackett  Burman  analyses  indicate  that  the  three  structural  parameters,  P0,  G+l,  and 
Mb,  are  significant  in  the  early  decomposition  phase,  since  the  significant  factor  effects  are  only 
present  in  the  early  stages  of  mass  loss  between  480  and  780  seconds.  These  results  are  depen¬ 
dent  on  the  selection  of  the  low  and  high  values  of  each  parameter.  The  parameters  with  a  small 
difference  between  the  high  and  low  values,  such  as  pressure,  are  shown  to  have  small  factor 
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Table  4.  Calculated  factor  effects  from  Plackett  Burman  Analysis 


Confidence  interval  (double-sided  probability  point  of  t-distribution):  80%  (1.33)  90%  (2.13)  95%  (2.78)  98%  (3.75)  99%  (4.60) 
Significant  factor  effects  based  on  90%  confidence  interval  using  four  degrees  of  freedom. 
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effects.  If  the  difference  between  the  high  and  low  values  were  greater,  the  factor  effect  may  have 
been  more  significant.  Specification  of  the  initial  bridge  population,  P0,  is  more  significant  than 
specification  of  the  initial  strong  bridge  population,  c0,  as  shown  in  Table  4.  The  lattice  coordina¬ 
tion  number,  a+1,  and  the  molecular  weight  of  bridges  are  also  significant  factors  that  describe 
the  chemical  structure  of  the  foam. 

The  Plackett  Burman  analysis  also  shows  that  all  of  the  steps  in  the  bond-breaking  scheme 
are  significant  except  for  Reaction  (7)  in  Table  1.  The  analysis  also  indicates  that  the  most  sensi¬ 
tive  pathway  is  the  formation  of  the  secondary  polymer.  At  the  expense  of  losing  generality,  the 
PUF  model  could  possibly  be  simplified  for  the  20  °C/min  TGA  simulation. 


7.2  GRID-INDEPENDENT  SOLUTIONS  USING  FINITE  ELEMENT  ANALYSIS 

Calculation  of  the  object  function  for  the  response  derivative  analysis  discussed  in  Section  7.3 
requires  grid-independent  solutions  using  a  finite  element  heat  conduction  code.  The  PUF 

decomposition  model  is  implemented  into  the  finite  element  heat  transfer  code,  COYOTE,1 
employing  a  user  subroutine.  COYOTE  is  a  two-  or  three-  dimensional,  finite  element,  massively 
parallelized  computer  program  designed  for  analysis  of  nonlinear  heat  conduction  problems.  In 
addition  to  solving  standard  thermal  diffusion  problems,  COYOTE  includes  the  effect  of  phase 
change,  condensed-phase  chemistry,  and  surface-to-surface  radiation.  Material  deletion  is 
included  through  the  use  of  the  finite  element  “death”  capability.  Elements  can  be  removed  to 
create  a  “dynamic”  radiation  enclosure.  Viewfactors  are  recalculated  whenever  an  element  dies. 
Gaps  at  material  interfaces  are  modeled  as  contact  surfaces.  Material  properties  can  be  tempera¬ 
ture  and/or  species  dependent  and  either  isotropic  or  orthotropic.  A  wide  variety  of  boundary 
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conditions  are  supported  in  COYOTE,  and  pre  and  postprocessing  file  formats  are  used  which 
permit  integration  with  existing  meshing  and  graphics  visualization  programs.  For  the  simula¬ 
tions  reported  in  the  current  work,  elements  are  dynamically  removed  from  the  simulation  when 
the  condensed  fraction  within  an  element  dropped  below  1%.  The  velocity  of  the  decomposition 
front,  referred  to  as  the  burn  front  in  this  report,  is  determined  from  the  time  of  element  death.  In 
The  decomposition  front  is  referred  to  as  a  “burn  front”  since  burning  is  defined  as  material 
destruction  by  heat.  The  foam  recession  rate  is  referred  to  in  this  report  as  the  foam  “burn  rate.” 

The  one-dimensional  simulation  of  burning  foam  was  performed  with  COYOTE  using  a 
single  column  of  elements  as  shown  on  the  left  side  of  Fig.  15.  Three  sides  of  the  column  of  ele¬ 
ments  were  assumed  to  be  insulated  and  the  fourth  side  was  exposed  to  various  boundary  condi- 

tions  such  as  radiation.  Experimentally  determined  thermal  conductivity  and  specific  heat  were 
used  in  the  one-dimensional  simulations.  Several  methods  can  be  used  to  determine  the  grid  sen¬ 
sitivity  of  decomposing  polyurethane  foam.  For  example,  the  size  of  the  element  can  be  reduced 
until  there  is  no  change  in  the  bum  front  velocity  as  shown  in  the  plot  of  Fig.  15.  Another  method 
would  be  to  specify  a  sufficient  number  of  elements  across  the  reaction  front  to  resolve  steep  gra¬ 
dients  in  temperature  or  species  profiles. 

The  grid  error  in  the  legend  of  Fig.  15  is  defined  as 

grid  error  =  100  x(Vi-V)/Vi,  (50) 

where  V  represents  the  burn  velocity  and  V,-  represents  the  grid-independent  bum  velocity.  The 

grid-independent  foam  recession  rate  in  Fig.  15  is  1.14  cm/min.  Figure  16  shows  the  grid  error 
for  a  one-dimensional  model  with  two  separate  constant  flux  boundary  conditions.  The  grid-inde- 
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.15  Dependence  of  foam  recession  rate  on  element  size  for  a  one-dimensional  model  with 
a  radiative  boundary  condition. 


Fig.  16  Grid  error  calculated  from  Eq.  (50)  for  a  one-dimensional  model  with  flux  boundary 
conditions. 
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pendent  foam  recession  rate,  V),  is  equal  to  0.64  cm/min  and  0.31  cm/min  for  the  1.25  cal/cm2-s 

and  0.6  cal/cm2-s  boundary  conditions,  respectively. 

In  Fig.  16,  the  predicted  velocity  using  2-mm  grids  is  0.54  cm/min,  giving  a  grid  error  of 

15%  for  the  1.25  cal/cm2-s  case.  Simulations  using  elements  smaller  than  2-mm  predict  faster 
burn  rates  than  0.54  cm/min.  The  faster  burn  rates  for  the  simulations  using  smaller  elements  is 
related  to  element  death.  For  example,  in  the  simulations  using  smaller  elements,  the  criteria  for 
element  removal,  condensed  fraction  being  less  than  0.01,  is  satisfied  earlier  than  the  criteria 
would  be  satisfied  for  a  simulation  using  larger  elements.  The  smaller  elements  are  removed 
from  the  computational  domain,  and  the  boundary  condition  is  then  applied  to  the  newly  exposed 
element.  For  an  element  size  of  0.5-mm,  the  grid  error  is  less  than  five  percent.  The  temperature 
gradient  across  each  element  is  also  an  indication  of  when  grid  independence  is  achieved.  For 
example,  grid  independence  in  the  present  study  was  achieved  when  the  temperature  gradient 
across  a  dying  element  was  approximately  50  °C  or  less.  Three-dimensional  calculations,  using 
0.05-mm  elements  are  currently  not  practical  without  adaptive  gridding.  However,  0.5-mm  ele¬ 
ments,  can  be  used  with  some  expense. 

7.3  RESPONSE  DERIVATIVE  ANALYSIS 

The  object  oriented  optimization  code,  DAKOTA,20  used  to  determine  the  PUF  model  activation 
energies,  can  also  be  used  to  determine  the  gradient  of  an  object  function  with  respect  to  model 
parameters.  This  response  derivative  can  be  used  to  quantitatively  determine  the  uncertainty  in 
predicted  results.  In  this  report,  two  response  derivative  analysis  were  performed  for  fast 
(0.97±0.14  cm/min)  and  slow  (0.26±0.05  cm/min)  burning  conditions.  For  the  response  deriva¬ 
tive  analysis,  the  object  function  was  chosen  to  be  the  steady-state  foam  recession  rate  for 
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one-dimensional  foam  decomposition,  as  discussed  in  Section  7.2.  The  TGA  simulations,  dis¬ 
cussed  in  Section  7.1,  did  not  require  solution  of  an  energy  equation,  since  the  temperature  history 
of  the  foam  samples  were  specified.  However,  the  regression  of  the  decomposition  front  depends 
on  the  thermophysical  properties  of  the  foam, 
as  well  as  the  decomposition  mechanism. 

Temperature  dependent  thermophysical  prop¬ 
erties  for  a  rigid  polyurethane  foam  are  given 
in  Table  5.  The  temperature  dependent  ther¬ 
mal  conductivity  and  the  heat  capacity  were 

measured23  with  a  differential  scanning  calo¬ 
rimeter  and  a  laser  flash  technique,  respec¬ 
tively. 

Two  different  rates  were  investigated: 

1)  a  fast  burn  rate  (0.97±0.14  cm/min) 
obtained  using  a  a  radiative  boundary  condi¬ 
tion  with  a  source  temperature  of  1,000°  C 
and  2)  a  slow  bum  rate  (0.26±0.05  cm/min) 
obtained  using  a  constant  flux  boundary  con- 

dition  of  0.5  cal/cm  -s.  All  the  simulations  in  this  section  were  made  using  0.5-mm  elements. 
The  gradient  of  the  object  function  with  respect  to  model  parameters,  dV/dC, was  used  to  deter¬ 
mine  the  uncertainty  in  foam  recession  rate  with  the  following  equation: 


Table  5.  Thermophysical  properties  of  foam 


Variable 

Symbol 

Description 

Nominal 

Value 

P  0 

Foam  density,  g/cm3 

0.353 

kf 

Thermal  conductivity,5 
cal/s-cm-K 

23°  C 

1.4x1  O'4 

50°  C 

1.5xl0“4 

100°  C 

1.6xl0"4 

150°  C 

1.8xl0“4 

200°  C 

2.0xl0"4 

250°  C 

2.2xl0“4 

cpf 

Heat  capacity,5  cal/g-K 

23°  C 

0.303 

50°  C 

0.324 

100°  C 

0.358 

150°  C 

0.440 

200°  C 

0.475 

250°  C 

0.526 

(AT)2 


I 

i  =  1 


jv  aCa2 

'  C  J  ' 


(51) 
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where  V,  Q,  AV,  and  A Q  represent  the  foam  recession  rate,  the  ith  model  parameter,  uncertainty  in 

foam  recession  rate,  and  the  uncertainty  in  the  ith  parameter,  respectively.  The  subscript  i  repre¬ 
sents  the  initial  density  of  the  foam,  the  conductivity  of  the  foam,  the  heat  capacity  of  the  foam, 
the  initial  temperature  of  the  foam,  the  activation  energies  of  the  nine  PUF  reactions,  the  decom¬ 
position  enthalpy,  and  the  emissivity  of  the  foam.  Since  the  thermal  conductivity  and  heat  capac¬ 
ity  of  the  foam  are  temperature  dependent,  the  sensitivity  parameter  for  these  variables  was 
modified  by  using  a  multiplying  factor  to  scale  the  temperature-dependent  property.  The  scale 
factor  for  both  the  thermal  conductivity  and  specific  heat  was  chosen  to  be  one.  A  similar 
approach  was  used  for  the  reaction  enthalpy. 

Table  6  gives  the  15  parameters  used  in  the  response  derivative  analysis  of  the  fast  bum 
rate  obtained  using  the  one-dimensional  burn  model  with  a  radiative  boundary  condition.  The 
nominal  parameter  values,  Q;  the  gradient  of  the  foam  recession  rate,  dV/cK^v  the  estimated 

parameter  uncertainty,  A  C,{,  and  the  individual  variable  uncertainty,  {y=rf  ;  are  also  given  in 

Table  6.  The  contributions  to  the  overall  uncertainty  is  shown  in  parenthesis  in  Table  6.  Vari¬ 
ables  with  a  contribution  to  the  overall  uncertainty  greater  than  5%  are  shaded.  The  most  sensi¬ 
tive  variables  affecting  foam  recession  rate  for  the  case  with  the  radiative  boundary  condition  are 
the  activation  energies  associated  with  the  Reactions  (2),  (4),  and  (5)  as  given  in  Table  1.  These 
are  the  same  kinetic  variables  that  were  determined  to  be  the  most  sensitive  kinetic  variables  in 
the  Plackett  Burman  analysis  discussed  in  Section  7. 1 .  A  one  percent  change  in  model  parameters 
resulted  in  a  14%  change  in  the  high  foam  recession  rate  case. 
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Table  6.  Uncertainty  Analysis  for  Radiation  Boundary  Condition 


i 

c  * 

dV/dQ 

AC,' 

rdV 

.m/  cJ 

,  cm2/min2 

Po 

0.353  g/cm3 

9.93X10'1 

0.0  K, 

0.0000123  (0.1%) 

kf 

1.0 

2.93x10° 

O.Olx^j 

0.0008601  (4.4%) 

cpf 

1.0 

1.90x10° 

O.OIXC; 

0.0003626(1.9%) 

To 

300  K 

-8.59xl0'3 

0.0  K, 

0.0006648  (3.4%) 

El 

53400  cal/mol 

-1.94x1 0'5 

O.Olx^j 

0.0001071  (0.6%) 

E2 

44580  cal/mol 

-9.90x1 0'5 

O.OlxC, 

0.0019478  (10.0%) 

e3 

42520  cal/mol 

7.33X10'6 

O.OlxC; 

0.0000097  (0.0%) 

e4 

45900  cal/mol 

-2.30x1  O'4 

O.OlxC, 

0.0111401  (57.3%) 

e5 

44600  cal/mol 

1.28X10'4 

O.Olx^j 

0.0032844  (16.9%) 

e6 

58060  cal/mol 

2.79xl0'5 

O.Olx^. 

0.0002632(1.4%) 

e7 

51850  cal/mol 

2.32X10'5 

O.Olx^j 

0.0001449  (0.7%) 

e8 

52830  cal/mol 

-2.31xl0“5 

O.OlxC, 

0.0001492  (0.8%) 

e9 

56660  cal/mol 

-3.10xl0'5 

O.Olx^- 

0.0003085  (1.6%) 

hrf 

-35  cal/cm3 

3. 49x1  O'2 

O.Olx^j 

0.0001492  (0.8%) 

£ 

0.8 

7.56xl0_1 

O.OlxC, 

0.0000366  (0.2%) 

£  =  0.01944  (100%) 

*  The  foam  recession  rate  for  a  radiation  boundary,  with  source  temperature  of 
1,000  °C,  is  0.97+0.14  cm/min. 


The  sensitivity  of  the  PUF  model  parameters  to  a  moderate  boundary  condition  of 

0.5  cal/cm2-s  was  also  investigated.  The  moderate  boundary  condition  resulted  in  a  thicker  reac¬ 
tion  zone  making  variables  such  as  thermal  conductivity  more  significant.  As  highlighted  in 
Table  7,  the  most  sensitive  variables  affecting  foam  recession  rate  for  the  constant  flux  calcula¬ 
tion  are  the  initial  density,  thermal  conductivity,  kinetics  effecting  the  formation  and  destruction 
of  the  secondary  polymer,  and  the  reaction  enthalpy.  A  one  percent  variation  in  the  input  factors 
resulted  in  a  18%  change  in  the  foam  recession  rate.  The  sensitivity  of  the  model  to  input  param- 
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Table  7.  Uncertainty  Analysis  for  Constant  Flux  Boundary  Condition 


i 

c,- 

dV/dQ 

AC,' 

rdV  2,  •  2 

c  f-rp — =-  >  cm  /mm 

dt*i  +'  > 

P  0 

0.353  g/cm3 

-4.51x10° 

O.OlxQ 

0.0002540  (11.1%) 

kf 

1.0 

1.09x10° 

O.OlxQ 

0.0001178  (5.1%) 

cpf 

1.0 

1.48X101 

O.OlxC,- 

0.0000022  (0.1%) 

T0 

300  K 

-3.09x10  3 

O.OlxQ 

0.0000862  (3.8%) 

Ei 

53400  cal/mol 

2.53xl0'6 

O.OlxQ 

0.0000018  (0.1%) 

e2 

44580  cal/mol 

-1.68x1  O'5 

O.OlxC,- 

0.0000561  (2.5%) 

Es 

42520  cal/mol 

-6.30x1 0"6 

O.OlxQ 

0.0000072  (0.3%) 

e4 

45900  cal/mol 

-5. 20x1  O'6 

O.OlxQ 

0.0000057  (0.2%) 

e5 

44600  cal/mol 

-4.16xl0"5 

O.OlxC,- 

0.0003446  (15.1%) 

e6 

58060  cal/mol 

-2.88xl0"5 

O.OlxQ 

0.0002800  (12.2%) 

Ey 

51850  cal/mol 

7.95xl0"6 

O.OlxQ 

0.0000170  (0.7%) 

e8 

52830  cal/mol 

-3.61X10'5 

O.OlxQ 

0.0003627  (15.9%) 

e9 

56660  cal/mol 

3.06xl0"5 

O.OlxQ 

0.0003012  (13.2%) 

hyf 

-35  cal/cm3 

6.07x1 0"2 

O.OlxC,- 

0.0004516  (19.7%) 

£  =  0.002288  (100%) 

*The  foam  recession  rate  for  a  constant  flux  of  0.5  cal/s-cm2  is  0.26+0.05  cm/min. 


eters  show  that  the  foam  decomposition  process  is  highly  nonlinear.  However,  since  the  foam 
recession  rate  gradients  contain  both  positive  and  negative  values,  the  18%  variation  in  foam 
recession  rate  is  a  conservative  estimate. 


8.  SCALE-UP  EXPERIMENT 

Chu  et  al.10  have  performed  component-scale  experiments  of  foam  decomposition  with 
well-characterized  boundary  conditions  using  a  heat  lamp  array  to  produce  fire-like  heat  fluxes. 
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In  Fig.  17,  the  two-dimensional  axisymmetric  geometry  of  the  foam  experiment  is  shown  with  an 
8.8-cm  diameter  by  15-cm  long  right  circular  cylinder  of  foam  encapsulating  a  3.8-cm  diameter 
by  6.4-cm  long  right  circular  cylinder  of  304  stainless  steel  (SS).  The  face  of  the  embedded  com¬ 
ponent  was  3.2-cm  from  the  heated  surface.  The  foam  cylinder  was  contained  in  a  sample  cup 
that  is  shown  in  Fig.  18.  The  cup  consisted  of  a  0.95-cm  thick  stainless  steel  bottom,  which  was 
press  fit  into  a  7.3-cm  long,  thin  wall  (0.5-mm)  stainless  steel  tube.  Six  6-mm  diameter  holes 
were  drilled  through  the  side  of  the  stainless  steel  tub,  near  the  cup  bottom,  to  vent  decomposition 
gases.  The  number  and  size  of  the  vents  were  chosen  to  keep  the  pressure  within  the  cup  near 
ambient  conditions  to  minimize  mass  transfer  effects.  Figure  19  shows  the  vent  holes  in  an  X-ray 

image  of  the  component-scale  experiment  after  15-min  exposure  to  an  incident  flux  of  25  W/cm2. 
The  bum  front  is  also  visible  in  Fig.  19. 

The  white  circles  in  Fig.  19  show  the  location  of  the  thermocouples  which  were  welded  to 
the  side  of  the  can  and  the  bottom  of  the  can.  The  location  of  the  thermocouples  on  the  can  are 
shown  in  more  detail  in  the  insert  of  Fig.  20.  Figure  20  also  gives  the  measured  temperature  at 
various  axial  can  locations  as  well  as  the  temperature  of  the  bottom  of  the  can  for  the  compo¬ 
nent-scale  experiment.  The  boundary  conditions  for  the  COYOTE  simulation  of  the  compo¬ 
nent-scale  experiment  were  set  to  the  specific  temperatures  measured  along  the  surface  of  the 
confining  cup.  The  cup  bottom  temperature  was  maintained  at  1,000°  C  with  the  thermocouple 
located  within  the  bottom  plate  as  shown  in  Fig.  20.  Figure  21  shows  the  measured  cup  wall  tem¬ 
peratures  and  bottom  cup  temperature  at  early  times.  In  Fig.  21,  temperatures  measured  with 

thermocouples  ©  and  ©  should  not  be  hotter  than  the  temperatures  measured  with 

thermocouple  o 
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Heat  lamp  array 


rxrxr^  Foam 


Fig.  17  Schematic  of  component- 


Slot  in  tube  for  bottom 
plate  thermocouple 

303  SS  bottom  plate 
press  fit  into  tube 

6-mm  diameter  weep  holes 
Cup  thermocouples 
321  SS  tube  (0.5-mm  wall) 


Fig.  18  Foam  cylindrical  holding  cup.  The  foam  (not  shown)  is  15-cm  long  and 
protrudes  from  the  cup. 
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Fig.  19  X-ray  image  of  component-scale  experiment  after  15-min  exposure  to  an 
incident  flux  of  25  W/cm2. 
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Figure  22  shows  the  axisymmetric  mesh  and  an  X-ray  of  the  foam  with  the  encapsulated 
304  SS  component.  A  3-mm  gap  was  visible  in  the  X-ray  taken  prior  to  the  experiment.  The 
3-mm  gap  was  modeled  as  a  radiation  enclosure.  The  element  dimensions  were  chosen  to  be 
5  mm  by  5  mm  to  maintain  at  least  three  to  four  elements  across  the  reaction  zone  as  shown  in 
Fig.  23.  The  mesh  in  Fig.  22  is  expected  to  give  a  grid  error,  as  defined  by  Eq.  (50),  less  than  5%. 
The  temperature  and  density  profiles  across  the  decomposition  front  are  shown  in  Fig.  23. C.  The 
thickness  of  the  front  estimated  from  Fig.  23. C  is  4  mm,  which  confirms  that  the  0.5-ram  element 
dimension  is  sufficiently  small.  Figure  23. D  also  shows  the  predicted  temperature  profile  and 
solid  fraction  profile  across  the  bum  front.  The  density  profile  shown  in  Fig.  23. C  has  significant 
scatter.  The  density  profile  was  estimated  from  the  pixel  density  in  the  X-ray  image  shown  in 
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Fig.  22  Axisymmetric  mesh  with  1 1,209  elements  and  X-ray  of  foam  in  cylindrical 
holding  cup  showing  3  mm  gap  between  face  of  component  and  foam.  The 
foam  is  15-cm  (6  inches)  long  and  protrudes  from  the  cup. 
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Fig.  23  A)  Measured  and  predicted  bum  front  shape  after  10-min  exposure  to  an  incident  flux 
~25  W/cm2,  B)  predicted  solid  fraction  contours,  C)  measured  temperature  and  density 
across  burn  front,  and  D)  predicted  temperature  and  solid  fraction  across  burn  front. 


Fig.  23. A.  Despite  this  uncertainty,  the  measured  density  profile  and  the  predicted  solid  fraction 
profile  across  the  burn  front  are  similar. 

Figure  24  shows  a  comparison  between  X-ray  images  and  calculated  temperature  profiles 

of  the  foam  with  the  cup  bottom  exposed  to  an  incident  flux  of  ~25  W/cm2.  The  boundary  condi¬ 
tions  were  set  to  the  specific  temperatures  measured  along  the  cup  surface  and  bottom  as  shown  in 
Fig.  20.  The  nominal  values  of  the  PUF  model  parameters  given  in  Table  2  and  the  tempera¬ 
ture-dependent  thermophysical  properties  given  in  Table  5  were  used  for  the  calculations  shown 
in  Fig.  24.  Resistance  across  the  foam  residue  on  the  face  of  the  component  was  also  accounted 
for  by  assuming  a  reduced  component  surface  emissivity  of  0.6.  The  reaction  enthalpy  was  cho- 
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Fig.  24  X-ray  images  (left)  and  calculated  temperature  profiles  (right)  of  foam  with  cup  bottom 
exposed  to  an  incident  flux  of  ~25  W/cm2. 
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sen  to  be  consistent  with  bond  energies  of  polyatomic  compounds.  In  Fig.  24,  the  decomposition 
front  is  horizontal  prior  to  reaching  the  component.  As  the  front  gets  closer  to  the  embedded 
component,  the  decomposition  front  curves  around  the  component.  The  exact  location  of  the 
decomposition  front  is  seen  as  a  density  variation  in  the  X-rays.  COYOTE’s  element  death 
option  was  used  to  remove  elements  when  the  condensed  fraction  within  individual  elements  was 
less  than  0.01.  The  shape  of  the  front  is  difficult  to  determine  near  the  wall  in  the  X-rays  because 
of  the  curvature  of  the  confinement.  Nevertheless,  the  calculated  and  measured  shape  of  the 
decomposition  front  appear  to  agree  very  well. 

Figure  25  shows  the  location  of  the  thermocouples  at  the  center  and  midradial  positions. 
Figures  26,  27,  and  28  show  a  comparison  between  predicted  and  measured  center,  midradial,  and 
component  temperatures,  respectively.  The  encapsulated  component  appears  to  heat  up  the  bulk 
of  the  foam  as  shown  by  the  temperature  profile  labeled  ®  in  Fig.  26.  In  previous  experiments 

without  an  encapsulated  component,10  all  temperature  measurements  remained  essentially  con¬ 
stant  until  the  decomposition  front  was  within  close  proximity.  With  the  encapsulated  compo¬ 
nent,  the  temperature  at  the  various  thermocouple  locations  is  shown  to  increase  gradually  until 
the  decomposition  front  is  within  close  proximity,  and  the  temperature  increase  is  more  rapid. 
This  behavior  is  shown  to  be  more  pronounced  in  the  experimental  temperature  profiles  plotted  in 
Figs.  26  and  27.  In  Fig.  26,  the  temperature  measured  by  thermocouple  ®  appears  to  be  100  °C 
higher  than  the  predicted  temperature  at  that  location.  This  difference  between  the  measurement 
and  calculation  is  thought  to  be  related  to  premature  exposure  of  the  thermocouple  sheath,  as 
shown  in  Fig.  29,  to  the  hot  radiating  can  and  subsequent  heat  transfer  by  conduction  to  the  ther¬ 
mocouple  tip.  Since  the  foam  burns  faster  along  the  wall  than  in  the  center  of  the  foam,  the  ther¬ 
mocouples  are  likely  to  become  partially  exposed.  The  exposed  thermocouples  can  become  hot 
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Fig.  25  A)  Center  thermocouple  locations  and  B)  midradial  thermocouple  locations. 


Fig.  26  Calculated  (solid  lines)  and  measured  (dashed  lines)  temperatures  at  center  ther¬ 
mocouple  locations  shown  in  Fig.  25. A. 
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Fig.  27  Calculated  (solid  lines)  and  measured  (dashed  lines)  temperatures  at  midradial 
thermocouple  locations  shown  in  Fig.  25. B. 


Fig.  28  Calculated  (lines)  and  measured  (symbols)  temperatures  of  the  encapsulated 
component  at  the  locations  indicated  in  the  inset. 
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A)  Top  view  of  exposed  thermocouple 


B)  Side  view  of  exposed  thermocouple 


Fig.  29  A)  Top  view  of  exposed  thermocouple  and  B)  side  view  of  exposed  thermocouple. 


65 


Scale-up  Experiment 


by  radiative  energy  exchange  with  the  hot  boundary.  Furthermore,  heat  transfer  between  the  hot 
tube  walls  and  the  thermocouple  leads  may  have  also  resulted  in  unaccounted  heat  transfer  to  the 
thermocouple  junctions. 

Figure  30  shows  the  predicted  transmitted  heat  flux  along  the  center  of  the  foam  at  various 
times.  The  maximum  transmitted  heat  flux  within  the  foam  at  each  time  plotted  in  Fig.  30  was 
located  at  the  burn  front.  Furthermore,  the  overall  maximum  transmitted  heat  flux  through  the 
foam  for  all  times  occurred  at  early  times  when  the  foam  surface  was  intimate  with  the  heated  sur¬ 
face.  As  the  burn  front  moved  away  from  the  heated  surface,  the  maximum  transmitted  heat  flux 
decreased.  The  decrease  in  maximum  transmitted  heat  flux  is  likely  due  to  view  factor  effect 
since  gas  opacity  and  convective  heat  transfer  was  not  included  in  the  model. 

Figure  31  shows  the  predicted  center  and  midradial  burn  front  velocities.  The  velocities 
were  estimated  using  the  time  of  element  death.  The  center  burn  front  velocity  was  slightly 
greater  than  the  midradial  burn  velocities,  resulting  in  a  curved  burn  front  shown  previously  in 
Fig.  24.  The  center  burn  front  velocity  profile  stops  at  about  7  minutes  when  the  burn  front 
reaches  the  embedded  component. 

The  2-D  axisymmetric  mesh  contained  1 1,209  elements  as  shown  in  Fig.  22.  The  simula¬ 
tion  required  11  days  and  16  hours  of  CPU  time  on  a  single  processor,  3  days  16  hours  of  CPU 
time  on  four  processors,  and  2  days  1  hour  on  eight  processors.  Ninety  percent  of  the  CPU  time 
was  attributed  to  the  view  factor  computation,  while  only  nine  percent  of  the  computation  was 
attributed  to  the  PUF  decomposition  model.  Adaptive  gridding  would  reduce  the  chemistry  com¬ 
putation  considerably  and  increase  simulation  accuracy,  but  the  view  factor  CPU  time  would 
likely  increase.  The  view  factor  calculation  would  dramatically  increase  if  the  adaptive  grid 
refinement  approaches  resolution  scales  defined  in  Section  7.2. 
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Fig.  31  Predicted  center  and  mid-radial  burn  front  velocity.  The  center  and  mid-radial 
locations  are  shown  in  the  inset. 
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9.  SUMMARY  AND  CONCLUSIONS 

The  PolyUrethane  Foam  (PUF)  decomposition  model  has  been  applied  to  thermal  decomposition 
of  rigid  polyurethane  foam  using  the  three-dimensional  finite  element  code  COYOTE.  The  PUF 
model  is  comprised  of  a  kinetic  mechanism  that  describes  bond  breaking,  a  lattice  statistics  model 
to  describe  the  mass  fraction  of  finite  polymer  fragments  isolated  from  the  macromolecule  by  bro¬ 
ken  bonds,  and  a  transport  model  to  describe  the  evaporation  of  the  finite  polymer  fragments  with 
high  vapor  pressures  into  the  gas  phase.  The  chemical  structure  of  the  rigid  polyurethane  foam 
was  determined  from  specific  synthesis  procedures,  and  the  PUF  model  parameters  were  obtained 
from  the  most  probable  structural  units  of  the  foam.  Kinetic  parameters,  for  the  PUF  bond-break¬ 
ing  mechanism,  were  obtained  by  using  an  object-oriented  optimization  code,  used  to  minimize 
the  difference  between  weight  loss  measured  in  a  small-scale  Thermal  Gravimetric  Analysis 

(TGA)  apparatus  using  nominally  5-mg  samples  from  six  isothermal  experiments10  and  model 

predictions.  The  temperature-dependent  thermal  conductivity  and  heat  capacity  were  measured23 
with  a  differential  scanning  calorimeter  and  a  laser  flash  technique,  respectively.  The  reaction 

enthalpy  was  obtained  from  a  differential  scanning  calorimeter.10  Three  nonisothermal  experi¬ 
ments,  not  used  to  estimate  kinetic  parameters,  were  simulated  to  validate  the  selection  of  the 
kinetic  parameters  for  the  PUF  model. 

Various  sensitivity  and  uncertainty  analysis  were  performed  using  the  PolyUrethane  Foam 
(PUF)  decomposition  model.  Although  exact  agreement  between  experimental  data  and  predic¬ 
tions  were  not  obtained  in  this  study,  uncertainty  in  both  experiment  and  model  have  been 
explored  and  quantified.  In  the  experiments,  a  few  of  the  thermocouples  used  to  define  the 
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boundary  conditions  were  shown  to  unexpectedly  increase  in  temperature.  Nonuniform  bum 
fronts  exposed  thermocouple  sheaths  to  the  hot  boundary.  Conduction  of  heat  along  the  thermo¬ 
couple  sheath  to  the  thermocouple  tip  may  have  caused  the  measured  thermocouple  to  read  higher 
temperatures  than  the  surrounding  foam.  Heat  transfer  between  the  hot  tube  walls  and  the  ther¬ 
mocouple  leads  may  have  also  resulted  in  unaccounted  heat  transfer  to  the  thermocouple  junc¬ 
tions.  The  uncertainty  in  the  exposed  thermocouples  are  on  the  order  of  100  °C.  Future 
experiments  will  use  more  thermocouples  to  better  defined  boundary  conditions. 

Two  separate  methods  were  used  to  determine  model  sensitivity  with  similar  results:  1)  a 
variance  analysis  and  2)  a  response  derivative  analysis.  The  variance  analysis,  based  on  experi¬ 
mental  design  techniques,  was  used  to  determine  variables  that  have  primary  effects  on  the  con¬ 
densed  fraction  predicted  using  a  20  °C/min  temperature  ramp.  The  primary  effect  variables 
included  the  initial  polymer  structure  and  the  kinetics  associated  with  the  formation  of  a  second¬ 
ary  polymer.  Two  response  derivative  analysis  were  performed  by  calculating  the  gradient  of  the 
foam  recession  rate  with  respect  to  model  parameters  for  a  fast  regression  (0.97±0.14  cm/min) 
obtained  using  fire-like  boundary  conditions  and  a  slow  regression  (0.26±0.05  cm/min)  obtained 
using  a  mild,  constant  flux  boundary  condition.  For  the  fast  regression,  the  most  sensitive 
kinetic  parameters  were  associated  with  the  formation  of  the  secondary  polymer  in  agreement 
with  the  Plackett  Burman  analysis.  For  the  slow  regression,  the  reaction  zone  thickened;  and  ther¬ 
mal  conductivity,  initial  density,  and  reaction  enthalpies  were  shown  to  be  sensitive  parameters  as 
well  as  the  activation  energies  associated  with  the  destruction  of  the  secondary  polymer.  The 
response  derivatives  were  significantly  different  for  the  fast  regression  case  than  the  slow  regres¬ 
sion  case.  Different  conditions  will  undoubtedly  give  different  sensitivity  results. 

A  two-dimensional  axisymmetric  simulation  of  an  unconfined,  component- scale  experi¬ 
ment  was  performed  using  an  1 1,209  element  mesh  with  the  PUF  model  implemented  into  COY- 
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OTE.1  The  simulation  of  decomposing  polyurethane  foam  around  an  encapsulated  component 
included  a  3-mm  gap  at  the  face  of  the  cylindrical  component.  Predictions  were  made  using  a 
grid-independent  solution  with  0.5-ram  by  0.5-mm  elements  and  compared  to  experimental  tem¬ 
perature  measurements.  Elements  were  removed  from  the  computational  domain  whenever  the 
condensed  fraction  was  less  than  0.01,  and  view  factors  were  recomputed.  Good  agreement  was 
obtained  between  temperature  measurements  and  predictions.  The  shape  of  the  regression  front 
measured  with  X-ray  imaging  was  also  in  agreement  with  the  predicted  shape  of  the  regression 
front.  Although  the  agreement  between  the  predicted  and  measured  temperatures  was  not  exact, 
the  temperatures  are  within  experimental  and  computational  uncertainty. 
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